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Abstract 

The  decision  making  model  (DMM)  [1,  2]  has  been  shown  to  generate 
phase  transitions,  to  be  topologically  complex  and  to  manifest  temporal 
complexity  through  size-dependent  random  fluctuations  in  the  switching 
times  between  the  two  critical  states  of  consensus.  These  properties  are 
entailed  by  the  fundamental  assumption  that  the  network  elements  in  the 
DMM  imperfectly  imitate  one  another.  The  process  of  subordination  es¬ 
tablishes  that  a  single  network  element  can  be  described  by  a  fractional 
master  equation  whose  analytic  solution  yields  the  observed  autocorrela¬ 
tion  function  obtained  by  numerical  integration  of  the  DMM  to  a  high 
degree  of  accuracy. 


1  Introduction 

Fractional  differential  equations  have  been  found  to  be  a  convenient  way  to 
describe  the  dynamics  of  complex  phenomena  that  are  characterized  by  multi¬ 
ple  scales,  that  is,  to  describe  the  time  evolution  of  fractal  processes  [3]  with 
no  single  scale  being  dominant.  In  spite  of  the  success  of  the  mathematical 
descriptions  of  such  processes  there  has  been  a  lack  of  identification  and  in¬ 
terpretation  of  mechanisms  that  entail  fractional  dynamic  equations  in  a  social 
context.  Herein  we  provide  an  explanation  of  one  source  of  a  fractional  dif¬ 
ferential  equation  that  describes  the  dynamics  of  a  complex  network  using  a 
fractional  master  equation. 

West  and  Turalska  [4]  adopted  imitation  as  a  basic  social/psychological 
mechanism  that  is  part  of  what  it  means  to  be  human.  Imitation  as  a  mech¬ 
anism  for  explaining  social  behavior  dates  back  to  the  turn  of  the  twentieth 
century  [5,  6,  7]  and  is  used  [4]  to  construct  a  complex  dynamic  social  network. 
In  this  social  network  the  hypothesis  that  imperfect  imitation  (termed  an  echo 
response)  is  a  fundamental  mode  of  human  behavior  is  made  and  implemented 
through  the  choice  of  interaction  coefficients  in  the  mathematics  of  the  decision 
making  model  (DMM)  [1,  2,  8].  The  resulting  dynamics  were  then  used  to  deter¬ 
mine  what  properties  emerged  as  a  consequence  of  the  echo  response  hypothesis. 
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It  was  shown  that  the  hypothesis  entails  behavior,  such  as  the  phenomenon  of 
consensus,  that  is  not  implicit  in  the  initial  statement  of  the  model.  In  other 
words  the  hypothesis  was  implemented  through  the  choice  of  the  coupling  co¬ 
efficients  in  the  DMM  and  the  resulting  critical  dynamics  was  used  to  explain 
such  collective  effects  as  herding,  flocking  and  schooling  as  well  as  manifesting 
stochastic  behavior.  The  DMM  is  sketched  out  in  Section  2. 

Stochastic  differential  equations  have  historically  been  derived  in  physical 
systems  using  Hamiltonian-based  arguments  in  which  the  system  of  interest  is 
coupled  to  an  inhnitely  large  environment.  These  derivations  produce  Langevin 
equations,  which  involve  coarse  graining  the  dynamic  description  so  that  the 
influence  of  the  environment  appears  in  the  equation  of  motion  as  a  linear  dissi¬ 
pation  and  random  fluctuations  that  are  interdependent  through  a  flluctuation- 
dissipation  relation  [9].  This  technique  basically  decomposes  the  system  into 
macroscopic  and  microscopic  time  scales  with  the  later  providing  the  random 
fluctuations  associated  with  thermal  effects.  More  recently  it  was  shown  that 
when  the  microscopic  dynamics  are  given  by  a  non-integrable  Hamiltonian  the 
separation  of  time  scales  no  longer  occurs  because  the  microscopic  time  scales 
diverge  to  macroscopic  size,  resulting  in  a  fractional  Langevin  equation  for  the 
system  dynamics  [10].  Herein  we  find  a  similar  result,  in  the  sense  that  the  dy¬ 
namics  of  the  smallest  component  of  the  complex  network,  the  individual,  can 
be  described  by  a  fractional  Langevin  equation  as  shown  in  Section  3.  However 
in  the  present  approach  no  coarse  graining  is  envoked. 

In  the  DMM  neither  the  single  individual  nor  the  dynamics  of  the  social  net¬ 
work  are  specified  by  Hamiltonians.  The  predicted  behavior  of  the  single  element 
dynamics  is  compared  with  the  numerical  results  from  the  DMM  implemented 
on  a  two-dimensional  lattice  in  Section  4.  The  influence  of  the  network  on  the 
dynamics  on  the  individual  is  measured  by  the  individual  autocorrelation  func¬ 
tion.  The  autocorrelation  function  is  shown  to  be  a  stretched  exponential  in  the 
subcritical  region,  to  be  a  modulated  stretched  exponential  in  the  supercritical 
domain  and  a  modulated  Mittag-LefBer  function  in  the  critical  region.  In  all 
three  regions  the  analytic  solutions  to  fractional  Langevin  equations  provide 
excellent  fits  to  the  numerical  calculations  of  the  complex  network. 

In  Section  5  some  conclusions  are  drawn. 


2  Decision  Making  Model  (DMM) 

The  DMM  of  a  complex  social  network  represents  the  dynamics  of  the  probabil¬ 
ity  for  an  individual  to  be  in  either  of  two  states,  yes  or  no,  up  or  down,  on  or 
off,  by  a  set  of  coupled  two-state  master  equations.  Turalska  et  al.  [2,  8]  showed 
that  this  set  of  equations  for  a  complex  network  has  the  same  structural  form 
as  that  of  the  individual  elements  of  the  network: 

d 

-^Pl  =  -912Pi  +  g2lP2,  (1) 

^,P2  = -g2iP2  + gi2Pi-  (2) 
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The  quantity  Pj  {t)  is  the  probability  of  the  network  being  in  the  state  j  =  1,2 
at  time  t  and  the  probability  is  normalized  at  all  times  such  that 

Pi{t)  +  P2it)  =  1.  (3) 

The  dynamics  are  determined  by  the  choice  of  the  functional  form  of  the  tran¬ 
sition  rates  in  Eqs.(l)  and  (2).  In  this  hrst  case  each  probability  is  influenced 
by  the  states  occupied  by  all  the  elements  of  the  network  as  determined  by  the 
transition  rates 


9ij{t)  =  go  exp  [K {T,j{t)  -  Ei(t)}]  ■,  i  ^  j  =  1,2  (4) 


where 


S.  it) 


Ns  it) 
N 


(5) 


N  denotes  the  total  number  of  elements  in  the  network  and  Ng  (t)  is  the  number 
of  elements  in  the  state  s  =  1,2  at  time  t.  The  parameter  K  is  the  control 
parameter  that  determines  the  strength  of  the  interactions  between  elements  of 
the  network.  It  is  evident  that  the  quantity  Eg  is  an  erratic  function  of  time, 
and  is  a  global  property,  obtained  from  an  observation  of  the  entire  network. 
We  define  the  stochastic  global  variable 


e(t)  =  Ei(t)-E2(t)  = 


Ni  it)  —  N2  it) 
N 


(6) 


whose  variability  is  characteristic  of  the  entire  network  of  echoes,  that  is,  the 
echoing  response  of  the  network  to  the  echoed  opinions  of  its  coupled  con¬ 
stituents. 

In  the  situation  where  the  number  of  nearest  neighbors  coupled  to  the  el¬ 
ement  of  interest  consists  of  all  the  other  individuals  in  the  network  we  have 
all-to-all  (ATA)  coupling.  Consider  the  ATA  coupling  case  and  assume  that 
the  total  number  of  elements  within  the  network  N  becomes  infinite.  In  the 
N  ^  00  case  the  fluctuation  frequencies  collapse  into  probabilities  according  to 
the  law  of  large  numbers  Eg  =  Ps-  In  physics  this  replacement  goes  by  the  name 
of  the  mean  field  approximation  in  which  case  the  transition  rates  in  the  ME 
are  written 

9^,  =  90  exp[-Kip.-p^)]  =  1,2.  (7) 

The  formal  manipulation  of  the  master  equation  even  in  this  simplified  case  is 
made  a  little  simpler  if  we  introduce  the  difference  in  the  probabilities 


n  =  Pi-P2-  (8) 

Subtracting  Eq.  (2)  from  Eq.  (1)  after  some  algebra  yields  the  highly  nonlinear 
rate  equation  for  the  difference  variable 

^'^  =  -i9i2+92i)'n  +  i92i-9i2)  (9) 
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where  the  nonlinearity  enters  through  the  transition  rate  dependence  on  the 
difference  variable 


912  =  50  exp  [-it'IT]  ;  521  =  5o  exp  [it'll]  (10) 

in  the  mean  field  approximation.  By  inserting  Eq.  (10)  into  Eq.  (9)  we  obtain 


dV 


(11) 


and  the  network  dynamics  are  determined  by  the  potential  function  V (IT),  which 
is  a  symmetric  double  well  potential  with  the  explicit  form 


y(n)  =  ^ 


K  +  l 

n  sinh  it'll - — —  cosh  it'll 

it 


(12) 


Note  that  the  social  network  is  not  described  by  a  Hamiltonian  and  yet  the 
global  dynamics  is  indeed  described  by  an  effective  Hamiltonian,  that  being  the 
double  well  potential. 

The  cooperative  behavior  of  the  infinitely  large  ATA  coupled  network  de¬ 
scribed  by  Eq.(ll)  is  that  of  an  overdamped  particle  hopping  from  one  potential 
minimum  to  the  other,  whose  position  is  H  within  the  potential  Eq.(12).  For 
A'  <  1,  half  of  the  nodes  are  in  the  state  1  and  half  are  in  the  state  2  because 
there  is  only  a  single  broad  minimum  in  the  potential.  At  the  critical  value 
of  the  control  parameter  K  —  Kc  =  1  a  bifurcation  occurs  and  the  potential 
develops  two  wells  separated  by  a  barrier  as  discussed  by  Turalska  et  al.  [2]. 
The  height  of  the  barrier  increases  with  the  value  of  the  control  parameter. 

It  is  interesting  that  at  the  critical  value  of  the  control  parameter  the  ATA 
version  of  the  DMM  undergoes  a  phase  transition.  Note  that  the  amplitude 
of  ^{t)  depends  on  the  value  of  the  control  parameter  K.  When  K  =  0,  all 
elements  in  the  network  are  independent  Poisson  processes;  thereby  an  average 
taken  at  any  moment  of  time  over  all  of  them  yields  zero.  Once  the  value  of  the 
coupling  becomes  nonzero,  K  >  0,  single  elements  are  less  and  less  independent, 
resulting  in  nonzero  averages.  The  quantity  Kc  is  the  critical  value  of  the  control 
parameter  K,  at  which  point  a  phase  transition  to  a  global  majority  state  occurs. 
In  numerical  calculations  we  use  the  time  average  =  (1C  (f)l)  as  a  measure  of 
this  global  majority.  More  precisely,  after  an  initial  10®  time  steps,  the  average 
is  taken  over  the  same  number  of  the  consecutive  time  steps  of  the  model.  In 
Figure  1  the  thin  line  indicates  the  ATA  phase  transition  as  measured  by  ^cg- 
The  other  phase  transitions  indicated  are  for  the  Ising  model  (dashed  line)  and 
the  DMM  model  on  a  two-dimensional  lattice  as  discussed  elsewhere  [8]. 

Real-world  networks  are  not  ATA  coupled  since  interactions  typically  have 
finite  range  and  elements  are  spatially  separated.  Moreover,  real-world  networks 
have  finite  numbers  of  elements.  It  is  therefore  useful  to  examine  how  strongly 
the  mean  field  solutions  are  violated  when  we  relax  these  constraints.  The 
stability  condition  can  be  violated  in  at  least  two  different  ways.  The  first  way 
is  by  reducing  the  number  of  elements  to  a  finite  value.  The  second  way  is 
by  restricting  the  number  of  links  so  the  network  no  longer  has  ATA  coupling. 
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Figure  1:  The  phase  diagram  for  the  global  variable  The  thin  solid  line  and 
the  dashed  line  are  the  theoretical  predictions  for  the  fully  connected  and  the 
two-dimensional  regular  network,  respectively.  In  both  cases  N  =  oo  and  the 
latter  case  is  the  Onsager  theoretical  prediction  [?]  for  a  two-dimensional  regular 
lattice.  The  dots  corresponds  to  the  global  states  observed  for  the  DMM  on  a 
two-dimensional  regular  lattice  {N  —  100  x  100  nodes)  and  go  —  0.01.  Periodic 
boundary  conditions  were  applied  in  the  DMM  calculations.  (From  [2]  with 
permission.) 


In  real-world  networks  both  sources  of  equilibrium  disruption  are  expected  to 
occur.  For  the  time  being  we  retain  the  ATA  coupling  within  the  network  and 
consider  the  number  of  elements  N  to  be  finite.  In  this  latter  case  we  can  no 
longer  make  the  mean  field  approximation  and  the  dynamic  picture  stemming 
from  the  above  master  equation  is  significantly  changed. 

If  the  number  of  elements  is  still  very  large,  but  finite,  we  consider  the 
mean-field  approximation  to  be  nearly  valid  and  replace  the  deterministic  equa¬ 
tion  Eq.(ll)  with  the  variable  11  replaced  by  the  global  variable  ^  to  obtain  a 
stochastic  differential  equation  [1,  2]: 


dHt) 

dt 


dV{0 

di 


+  £  W 


(13) 
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and  the  additive  fluctuations  e  (t)  have  amplitudes  that  are  computationally 
determined  to  be  on  the  order  of  l/iV  .  Here  again  we  see  the  influence  of  the 
imperfect  imitations  that  echo  between  individuals  in  the  network. 

Note  that  the  double- well  potential  in  the  mean  field  approximation  persists 
in  the  present  description.  The  random  fluctuations  induce  transitions  between 
the  two  states  of  the  potential  well.  Consequently,  for  a  network  with  a  finite 
but  large  number  of  elements  the  phase  synchronization  of  Eq.(ll)  is  not  stable 
and  the  stochastic  differential  equation  represents  the  dynamics  of  the  network 
that  must  be  solved.  Furthermore  the  fluctuations  can  drive  the  particle  from 
one  well  of  the  potential  to  the  other  when  its  amplitude  is  sufficient  to  traverse 
the  barrier  separating  the  wells.  However,  here  the  fluctuations  arise  from  the 
finite  number  of  elements  in  the  network  rather  than  from  non-existent  thermal 
excitations  and  are  consequently  suppressed  as  the  network  size  increases. 


000  2.50x10‘  5.00x10'  7.50x10*  1.00x10® 


1.0 
0.5 
0.0 
-0.5 
-1.0 

000  2.50x10*  5.00x10*  7.50x10*  1.00x10* 


'  1  '  1  '  1  '  ; 

— 

- ' - 1 - ' - 1 — 

- 1 - ' - 

Time,  t  [arb.  units] 


Figure  2:  The  fluctuation  of  the  mean  field-average  phase  as  a  function  of  time. 
For  a  system  of  =  500  elements  {top),  N  —  1500  elements  {middle),  and 
N  —  2500  elements  {bottom).  The  coupling  constant  is  if  =  1.05  and  g  =  0.01 
in  all  three  cases.  (From  [2]  with  permission.) 

Although  Fq.  (13)  is  written  in  the  continuous  time  representation,  in  prac¬ 
tice  the  numerical  calculations  of  DMM  correspond  to  the  adoption  of  a  finite 
integration  time  step  At  =  1.  Note  that  the  stochastic  rate  equation  Fq.(13) 
replaces  Fq.  (11)  in  the  case  of  a  finite  N,  and  that  Fq.  (11)  is  recovered  in  the 
ideal  case  N  =  oo.  Consider  the  ATA  coupling  condition  with  a  finite  number 


6 


of  elements  by  numerically  integrating  the  master  equation  for  each  element  in 
the  network  and  then  calculating  the  number  of  elements  in  each  of  the  two 
states.  In  Figure  2  the  fluctuating  global  variable  ^{t)  is  depicted  as  a  function 
of  time,  under  differing  conditions.  Notice  that  with  increasing  N  the  fluctua¬ 
tion  ^{t)  become  more  distinctly  dichotomous-like,  with  an  increasingly  sharp 
transition  from  the  ’up’  to  the  ’down’  state.  This  pattern  corresponds  to  the 
entire  network  keeping  a  decision  for  a  longer  and  longer  time  as  the  size  of  the 
network  increases.  The  condition  of  a  decision  lasting  forever  is  reached  in  the 
ideal  case  N  —  oo.  The  global  variable  fluctuates  between  the  two  minima  of 
the  double- well  potential  as  described  by  Eq.(13)  for  K  =  1.05  >  and  three 
values  of  the  size  of  the  network  corresponding  to  an  ever  increasing  influence 
of  the  echo  network.  The  single  element  follows  the  fluctuations  of  the  global 
variable,  switching  back  and  forth  from  the  condition  where  the  -fl  state  is 
preferred  statistically  to  that  where  the  —1  sate  is  preferred  statistically.  The 
complete  properties  of  the  ATA  DMM  are  explored  by  Turaska  et  al.  [2,  8]. 

3  Subordination  and  Fractional  Dynamics 

Note  that  if  attention  is  concentrated  on  a  single  network  element  when  the 
network  is  in  a  consensus  state  that  individual  would  still  appear  to  make  tran¬ 
sitions  according  to  an  exponential  distribution  as  suggested  by  Figure  3.  That 
is  the  situation  on  a  large  scale  view  of  the  calculations  on  a  two-dimensional 
lattice  with  nearest  neighbor  interactions  and  a  critical  control  parameter  value 
of  Kc  =  1.7.  The  strict  exponential  is  indicated  by  the  black  dotted  curve.  The 
single  particle  survival  probabilities  do  not  look  too  much  different,  the  light  gray 
dashed  curve  with  the  subcritical  value  AT  =  1.5  <  is  very  close  to  the  ex¬ 
ponential.  The  remaining  single  particle  curves,  whether  critical  K  =  Kc  =  1.7 
or  supercritical  K  >  Kc  appear  to  be  exponential  on  this  graph.  The  difference 
in  the  behavior  of  the  individual  from  that  in  the  non-interacting  state  would 
be  that  she  tends  to  be  more  reluctant  to  change  her  mind. 

The  deviation  of  the  individual  survival  probability  from  the  exponential 
form  appears  to  be  modest  when  compared  with  the  dramatically  greater  devi¬ 
ation  of  the  survival  probability  of  the  global  variable  from  the  exponential.  The 
average  network  behavior  differs  markedly  as  the  control  parameter  is  increased 
from  the  subcritical  through  the  supercritical  regions.  However  the  influence  of 
the  global  variable  on  the  behavior  of  the  individual  does  not  appear  to  induce 
a  significant  change.  For  the  individual  the  change  is  a  subtle  yet  profound 
difference  and  is  a  direct  result  of  the  hypothesized  imitation  mechanism.  So 
if  the  individual  survival  probability  is  not  exponential,  what  is  it?  To  answer 
this  question  we  turn  our  attention  to  describing  an  alternate  construction  of 
the  individual  element  dynamics. 

The  master  equation  for  a  single  isolated  element  is  according  to  the  numer¬ 
ical  simulation  given  by 

(j)  (nAr)  —  1]  At)  =  — {[n  —  1]  Ar)  At,  (14) 
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whose  discrete  solution  is 


^(n)  =  (1  -  goAT)”(/>(0) .  (15) 

Here  (p  (n)  is  the  difference  variable  for  a  single  individual  chosen  from  the 
network  at  random  and  as  n  ^  oo  and  At  ^  0  such  that  clock  time  is  t  =  nAr 
we  have  the  apparently  trivial  result 

(j){t)  =  .  (16) 


Figure  3:  The  survival  probability  dl  (t)  of  the  global  variable  (solid  lines)  is 
compared  with  the  transitions  between  two  stable  states  for  a  single  unit  Sj 
(dashed  lines) .  Simulations  were  performed  on  a  lattice  of  size  A  =  50  x  50  for 
go  —  0.01  and  increasing  values  of  the  control  parameter  K.  The  critical  value  of 
the  control  parameter  is  =  1.72.  Dotted  line  corresponds  to  the  exponential 
decay  with  rate  go-  The  slope  of  the  inverse  power  law  region  is  /r  —  1  =  0.50. 

However  when  the  individual  is  part  of  the  network  the  exponential  is  not 
what  is  observed.  The  straightforward  interpretation  of  the  index  n  given  above 
is  no  longer  correct.  Instead  there  exists  a  stochastic  connection  between  what 
is  now  the  discrete  operational  time  n  and  the  chronological  or  clock  time  t, 


one  that  is  referred  to  as  subordination  in  the  mathematics  literature  [11].  We 
adopt  the  subordination  interpretation  here  and  assume  that  the  time  interval 
t  lies  between  the  values  (n  —  1)At  <  t  <  nAr.  Consequently  the  time  t  is 
derived  from  a  waiting-time  pdf  given  by  that  of  the  network  as  a  whole: 


V’(f)  = 


(/z-  i)r^^i 
{T  +  tf 


(17) 


where  the  survival  probability  as  determined  from  Figure  3  to  be  [1,  2] 


•^{1)  =  jf) 

t 


rjn  \  /J,~l 

T^t) 


(18) 


An  event  occurring  at  time  n  is  an  abrupt  change  from  (f){n—l)  to  the  state 
(j)  (n) .  At  a  generic  time  t  we  can  write  the  average  probability  difference  as  the 
generalized  master  equation  (GME) 


{t’)^ 


(19) 


where  fin)  is  given  by  Eq.(15).  Here  we  see  that  the  GME  replaces  the  two- 
state  master  equation  of  the  DMM. 

Note  that  dtip^  (t)  is  the  probability  that  n  events  have  occurred,  the  last 
one  occurring  in  the  time  interval  {t,t  +  dt).  The  function  d/  (t)  denotes  the 
probability  that  no  event  occurs  up  to  time  t  and  is  given  by  Eq.(18).  The 
occurrence  of  an  event  corresponds  to  activating  a  decision  with  (1  —  goAr), 
so  that  activating  n  such  events  transforms  the  initial  condition  f  (0)  into  the 
product  (1  —  goAr)”  0  (0). This  form  of  the  equation  is  kept  from  time  t' ,  at 
which  time  the  last  event  occurs,  up  to  time  t,  the  time  interval  t  —  t'  being 
characterized  by  no  event  occurring.  Of  course,  the  expression  Eq.(19)  takes 
into  account  that  the  number  of  possible  events  may  range  from  the  no-event 
case  to  that  of  infinitely  many  events.  The  conditions  necessary  for  this  result 
to  occur  are  discussed  by  Svenkenson  et  al.  [12]. 

It  is  useful  to  introduce  Laplace  variables  in  our  discussion.  The  Laplace 
transform  of  a  function  f{t)  is  denoted 


7(s)  =  Je  ’'*f{t)dt.  (20) 

0 

Consequently  using  the  convolution  theorem  Eq.(19)  can  be  expressed  in  terms 
of  Laplace  transformed  quantities 

OO 

(f  (s))  =  (s)  ^  (s)  (1  -  ffoAr)”  4>  (0) .  (21) 

n—0 
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We  assume  the  intervals  between  successive  transitions  are  independent  of  one 
another  so  that  the  waiting  time  pdf  for  n  transitions  is  the  product  of  n  single 
transition  pdf’s: 

^  in 


V’n  (S) 


fis) 


(22) 


and  the  Laplace  transform  of  the  survival  probability  is 


$(s) 


(23) 


Inserting  these  last  two  expressions  into  Eq.(21)  and  evaluating  the  sum  yields 


- - - 

s  +  5oAr$  (s) 


(24) 


whose  inverse  Laplace  transform  yields  the  GME: 


=  -5oAr  t')  if  it'))  dt'. 


(25) 


3.1  Fractional  Langevin  Equation 


The  function  $  (t)  is  a  memory  kernel  containing  the  information  on  how  the 
other  elements  in  the  network  influence  the  dynamics  of  the  individual  element 
selected.  The  Laplace  transform  of  the  memory  kernel  is 


$(s) 


sip  (s) 

1  -  ^  (s) 


(26) 


and  a  complete  discussion  of  its  properties  is  now  given  in  textbooks  [3]. 

Previous  analysis  has  shown  that  the  waiting-time  pdf  is  an  inverse  power- 
law  distribution.  The  asymptotic  behavior  of  the  GME  is  determined  by  con¬ 
sidering  the  waiting-time  pdf  given  by  Eq.(17)  as  s  — >  0  : 


ip{s)  « 

so  that  Eq.(24)  reduces  to 


1  -  r(2  -  p) 


1 

S  + 


cPiO) 


and  the  rate  parameter  has  the  value 

50  At 


= 


r  (2  -  /z) 


1  <  A  <  2. 


(27) 


(28) 


(29) 


We  now  assume  that  the  exact  equation  for  the  single  person  dynamics  has  both 
an  average  and  a  fluctuating  part  just  as  we  found  in  the  mean  field  treatment  of 
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the  double  well  potential  [1,  8].  Consequently  in  terms  of  the  Laplace  variables 
we  have  the  stochastic  equation 


$(s) 


-e  (s) , 


which  has  the  inverse  Laplace  transform 


(30) 


This  is  a  stochastic  fractional  master  equation  or  fractional  Langevin  equation 
(FLE)  in  which  the  stochastic  properties  of  e  (t)  are  determined  by  the  fluctua¬ 
tions  of  the  complex  social  network.  The  fractional  derivative  in  this  equation 
is  of  the  Riemann-Liouville  form  [?] 


Dnm]^ 


1  V  f{t')dt' 

T  {m  —  q)  J  (t  — 

0 


(32) 


due  to  the  global  dynamics  of  the  network  q  =  n  —  1  >  0,  with  the  requirement 
that  m  is  the  smallest  integer  such  that  m—  1  <  q  <  m  so  that  m  =  1  in 
Eq.(32). 

The  solution  to  the  FLE  is  given  by 


<P{t) 


</)(0)i?^_i(-  (At)''”') 

t 

+Jit-  t'T'^  (-  (A  [t  -  i'])"'’')  £  {t')  dt'.  (33) 

0 


where  the  homogeneous  solution  to  the  fractional  equation  is  the  Mittag-Lefiler 
function  (MLF) 

k=0 

and  the  kernel  of  the  integral  is  in  terms  of  the  second  order  MLE  taken  up  in 
the  next  section. 

It  is  useful  to  replace  Eq.(30)  with  the  Laplace  transform  of  the  single¬ 
element  autocorrelation  function 


E 

(^(s))0(O) 

1 

E 

0(0)^' 

s  +  A'' 

(35) 


where  i?[-]  is  an  average  over  an  ensemble  of  initial  states.  Again  assuming 
the  existence  of  fluctuations  and  rearranging  the  terms  in  Eq.(35),  the  inverse 
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Laplace  transform  [?]  of  the  resulting  equation  yields  the  fractional  differential 
equation  for  the  autocorrelation  function 


Dr 


[cm 


r(2-M) 


=  + 


E  <j){0Y 


(36) 


It  is  clear  that  the  solution  to  the  homogeneous  fractional  equation  is  again  the 
MLF,  but  the  behavior  of  the  individual  autocorrelation  function  also  depends 
on  the  the  properties  of  the  random  fluctuations  as  we  show  in  the  next  section. 


4  Fractional  Dynamics  and  Complex  Networks 

The  solution  to  the  GME  is  obtained  from  the  exact  Laplace  transform  equation 
Eq.(24).  However  it  is  notoriously  difficult  to  obtain  analytic  expressions  by 
direct  inversion  of  the  resulting  eqnations  due  to  the  complexity  of  the  exact 
form  of  the  Laplace  transform  memory  kernel.  Consequently,  the  strategy  is  to 
consider  the  asymptotic  forms  of  the  solution,  which  we  do  by  examining  the 
solutons  to  the  FLE.  We  find  that  the  properties  of  the  fluctuations  change  as 
the  control  parameter  is  varied  from  the  subcritical,  critical  and  supercritical 
regions. 

4.1  Subcritical  Solution 

The  analytic  solution  to  Eq.(31)  with  the  noise  set  to  zero  is  obtained  elsewhere 
[?]  by  inverse  Laplace  transforming  Eq.(28).  The  solution  is  the  MLF,  as  we 
said, 

(37) 

The  autocorrelation  function  for  the  single  individual,  the  solution  to  the  frac¬ 
tional  differential  equation  Eq.(36),  is  therefore  also  given  by  the  MLF 

C{t)  =  (-  .  (38) 

From  the  series  form  of  the  MLF  it  is  evident  that  for  /i  =  2  the  autocor¬ 
relation  function  would  be  an  exponential.  Consequently  the  influence  of  the 
network  on  the  behavior  of  the  individual  in  this  case  would  be  essentially  that 
of  uncorrelated  random  noise  and  therefore  would  not  qualitatively  change  from 
the  Poisson  natnre  of  an  isolated  individual.  However  this  is  not  the  case  for 
other  values  of  the  inverse  power-law  index.  The  dynamic  behavior  of  the  net¬ 
work  results  in  a  stretched  exponential  autocorrelation  for  the  dynamics  of  the 
individual 


limC'(t)  =  1  — 

t-tO 


r(M) 


exp 


[X^\ 

r(t^)  )' 


(39) 
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Note  that  the  early  time  behavior  of  the  MLF  is  the  indicated  stretched  expo¬ 
nential. 

The  ten  individual  elements  used  in  the  evaluation  of  the  single  autocorre¬ 
lation  functions  were  selected  at  random  from  a  DMM  network  calculation  on 
a  two-dimensional  lattice  of  size  100  x  100  with  nearest  neighbor  interactions. 
It  is  clear  from  the  autocorrelation  functions  depicted  in  Figure  4  that  in  the 
subcritical  regime  the  data  are  not  indicative  of  exponential  decay.  For  the 
subcritical  case  oi  K  —  1.50  the  stretched  exponential  function  in  Eq.(39)  gives 
a  remarkably  good  fit  to  the  autocorrelation  of  the  time  series  data.  The  MLF 
index  is  seen  to  fall  in  the  interval  1.43  <  p.  <  1.55,  the  rate  of  decay  of  the 
stretched  exponential  is  in  the  range  3.4  x  10^^  <  A  <  4.8  x  10~^,  and  the 
quality  of  fits  all  lie  in  the  interval  0.97  <r^<  0.99.  The  fitting  of  the  analytic 
autocorrelation  function  at  early  time  to  the  numerically  generates  curves  using 
the  DMM  is  certainly  very  good  in  the  subcritical  domain. 


Figure  4:  Ten  realizations  of  the  single  element  autocorrelation  function,  each 
selected  randomly  from  the  elements  on  a  two-dimensional  lattice,  under  DMM 
dynamics  with  N  —  100  x  100,  —  0.01  and  K  —  1.50  are  plotted  with  thin  gray 

lines.  Thick  dashed  lines  denote  the  fit  with  a  stretched  exponential  function. 
The  hts  to  the  most  upper  and  lower  single  element  autocorrelation  functions 
are  shown. 
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4.2  Supercritical  Solutions 

In  the  supercritical  region  of  the  control  parameter  the  single  elements  used  in 
the  evaluation  of  the  autocorrelation  function  were  selected  at  random  from  the 
same  DMM  network  calculations  as  those  above.  The  calculations  were  done 
on  a  100  X  100  two-dimensional  lattice  with  nearest  neighbor  interactions.  It  is 
evident  from  the  autocorrelation  functions  depicted  in  Figure  5  that  there  exists 
a  wide  spread  in  either  the  kind  of  distributions  that  characterize  the  network 
dynamics  in  this  region  or  there  is  an  equally  wide  spread  in  the  parameter 
values  required  to  fit  the  data.  Let  us  consider  the  last  alternative  first. 

We  temporarily  assume  that  the  network  dynamics  in  the  supercritical  region 
is  similar  to  that  in  the  subcritical  region.  We  mean  by  this  that  the  stretched 
exponential  can  be  used  to  describe  the  autocorrelation  function  in  both  domains 
and  therefore  each  of  the  data  curves  in  Figure  5  can  be  fit  using  Eq.(39)  albeit 
with  a  wide  distribution  of  rates  and  exponents.  The  quality  of  the  fit  is  > 
0.84  for  each  of  the  curves  in  Figure  5,  which  is  remarkably  good  given  the 
vertical  spread  in  the  distal  values  of  the  autocorrelation  function  curves.  In 
many  fits  to  social  data  the  would  be  sufficient  to  draw  a  positive  conclusion. 
However  here  we  investigate  the  influence  of  the  noise  produced  by  the  network 
dynamics. 

We  assume  that  the  random  fluctuations  £  (t)  in  the  FLE  have  Gaussian 
statistics  and  are  exponentially  correlated  with  the  initial  condition  such  that 


E[<^{Q)e{t)]^2E 

In  addition  we  insert  the  identity 


(40) 


dE., 


_ V _ L  —  _\M“i 

dt 


(-  (At) 

into  Eq.(33)  to  obtain  for  the  autocorrelation  function 

(-  [A  (t-t')]""') 


fl-1 


C'(t)=E,_i(-  (At)''”^)-^ 


dt' 


e-^*'dt',  (41) 


Integrating  Eq.(41)  by  parts  and  using  the  method  of  steepest  decent  on  the 
remaining  integral  yields  for  the  early  time  solution 


limC'(t) 

t^o 


Dj 


A''- 


T+  1- 


D 


=  -7* 


\ 

1  exp 

r(/r)  _ 

(5^e^(‘‘)e-^‘erf  [y^2|r' 


A^ 


(42) 
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Figure  5:  The  autocorrelation  functions  (thin  gray  lines)  are  calculated  for  each 
of  100  time  series  obtained  from  a  DMM  calculation  on  a  two-dimension  lattice 
for  a  control  parameter  greater  than  the  critical  value  K  —  1.9  >  K^-  There  are 
ten  time  series  for  each  of  ten  individuals  chosen  at  random  from  the  network. 
Thick  dashed  lines  denote  the  fit  with  Eq.  (44) .  The  most  upper  and  lower  single 
node  autocorrelation  function  has  been  fitted. 


and  the  new  functions  are 


Fit*) 


7 


r(M) 


1) 

We  approximate  the  early  time  solution  to  be 


;  F"it*)  =  i2-n) 


7r  i^l) 


limC'(t)  «  — ? 

i-»0  ^  ^ 


\ 

1  exp 

r(Ai)  _ 

(43) 

(44) 


where  the  last  term  in  Eq.(42)  has  been  dropped.  We  subsequently  use  the 
fitted  parameter  values  to  calculate  the  size  of  the  neglected  term,  which  is 
determined  to  be  orders  of  magnitude  smaller  than  the  terms  retained  and 
therefore  its  neglect  is  justified. 
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Equation  (44)  is  used  to  fit  all  ten  curves  in  Figure  5  using  the  exponent 
/r  =  1.5  and  dissipation  rate  A  =  0.02.  The  quality  of  the  fit  to  all  the  curves 
on  the  graph  for  the  empirical  strength  of  the  fluctuations  D  and  the  decay 
rate  of  the  fluctuations  7  for  the  individuals  depicted  in  Figure  ??  fall  in  the 
range  0.96  <  <  0.99.  For  weakly  decaying  fluctuations  the  individual  au¬ 

tocorrelation  is  high  as  is  the  strength  of  the  fluctuations;  this  is  the  region 
of  exponentially  truncated  stretched  exponential  autocorrelation.  For  rapidly 
decaying  fluctuations  the  individual  autocorrelation  decays  as  a  stretched  expo¬ 
nential.  The  fitting  of  the  approximate  analytic  solution  to  the  autocorrelation 
function  to  the  numerically  generated  curves  using  the  DMM  is  certainly  very 
good  in  the  supercritical  domain. 


4.3  Critical  Solutions 


As  the  critical  point  is  approached  from  the  subcritical  region  the  autocorre¬ 
lation  function  changes  as  would  be  expected  due  to  the  formation  of  clusters 
as  the  network  undergoes  a  phase  transition.  The  plunging  stretched  exponen¬ 
tial  that  was  observed  in  the  subcritical  region  as  seen  in  Figure  4  is  replaced 
with  a  more  gently  decreasing  correlation.  Figure  6  depicts  another  ten  ran¬ 
domly  selected  single  element  trajectories  from  DMM  calculations  on  the  two- 
dimensional  lattice.  In  this  case  the  control  parameter  has  the  critical  value 
1.70  and  the  autocorrelation  function  is  calculated  from  the  single  particle  tra¬ 
jectories.  It  is  evident  by  comparing  these  data  with  the  data  curves  in  Figure  4 
that  the  autocorrelation  does  not  decrease  as  quickly  in  time  and  there  is  more 
variability  among  the  ten  slopes  at  late  time. 

The  question  is  whether  the  MLF  solution  to  the  FLE  can  be  used  to  fit  these 
data  as  well  as  it  did  the  suberitical  data  or  whether  we  require  an  additional 
mechanism  such  as  the  fluctuations  in  the  supercritical  region.  We  find  that  at 
early  times  the  stretched  exponential  approximation  to  the  MLF  fits  the  data 
very  well  as  shown  in  Figure  6.  At  late  times  correlated  fluctuations  generated 
by  the  remainder  of  the  network  must  be  included. 

The  data  curves  in  Figure  6  are  extremely  well  fit  using  the  stretched  ex¬ 
ponential  funetion  at  early  times,  however,  the  inverse  power  law  at  late  times 
using  the  same  MLF  is  not  a  good  fit.  The  fits  by  separate  asymptotic  solutions 
to  the  data  in  the  critical  region  are  excellent  each  with  >  0.99.  However 
there  is  a  problem  and  that  is  the  low  value  of  the  inverse  power-law  index, 
it  is  1.22  rather  than  the  1.5  found  in  the  numerical  calculation  of  the  sur¬ 
vival  probability  for  the  network  and  used  in  the  subordination  process.  This 
change  of  inverse  power-law  index  suggests  that  the  dynamic  equation  for  the 
autocorrelation  function  is  given  by  Eq.(30)  but  the  noise  is  no  longer  given 
by  exponentially  decaying  fluctuations.  The  strength  of  the  fluctuations  at  the 
critical  points  are  assumed  to  decay  as  an  inverse  power  law: 


E[cl,{0)e{t)]=E  <P{0) 


T  +  t 


(45) 
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Figure  6:  Ten  realizations  of  the  single  element  autocorrelation  function  selected 
randomly  from  the  two-dimensional  lattice  with  N  =  100  x  100,  go  —  0.01  and 
K  —  1.70  «  Kc-  Thick  dashed  lines  denote  the  fit  with  an  inverse  power  law 
function.  The  most  upper  and  lower  single  node  autocorrelation  function  has 
been  fitted..  The  data  are  fit  to  a  stretched  exponential  at  early  time  (dashed 
light  gray  line)  with  Aq  «  0.01  and  to  an  IPL  at  late  times  (dashed  dark  gray 
line)  with  index  /r  —  1  «  0.22.  The  quality  of  fit  >  0.99  for  all  ten  curves  in 
both  regions. 


so  that  the  Laplace  transform  of  the  autocorrelation  function  is 
1  D 


C{s)  = 


'  +  s  +  A''”^s2-a‘  (sT) 


(46) 


The  time  asymptotic  autocorrelation  function  is  obtained  by  considering  the 
s  — >  0  limit  and  noting  1  >  2  —  /r  as  well  as  —  1  <  0, 


limC  (s)  =  - — 

s^O  S  + 


1  +  r  {1  -  v)  D 


r(i-^)Dr-\,^,_3 

AM-1 
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and  consequently  yields  the  inverse  power  law  with  a  modihed  index 


lim  C(t)  =  - 

t— too  p 


T{l-iy)D 


f,-v){XTY 


T, 


fi+v—2 


(47) 


From  the  hts  to  the  curves  in  Figure  6  we  determine  that  the  average  fluctuation 
relaxation  rate  to  be  «  0.72  with  a  quality  of  fit  >  0.99  for  each  of  the 
ten  curves  considered. The  relaxation  properties  of  the  fluctuations  are  therefore 
seen  to  dominate  the  autocorrelation  of  the  individual  elements  of  the  network 
at  late  times  and  the  analytic  fit  in  the  two  time  regimes  are  uniformly  excellent. 


5  Conclusion 

The  numerical  solution  of  the  DMM  on  a  100  x  100  lattice  with  elements  at 
each  of  the  nodes  and  with  nearest  neighbor  interactions  gives  rise  to  an  inverse 
power-law  survival  probability  for  the  global  variable  introduced  in  Section  2. 
Using  the  notion  of  a  subordination  time,  that  being  the  time  experienced  by 
an  individual,  with  the  influence  of  the  network  entering  into  the  individual’s 
dynamics  through  the  distribution  of  critical  events,  the  dynamics  of  an  in¬ 
dividual  is  determined  by  a  fractional  master  equation.  The  explicit  form  of 
the  FLE  depends  on  whether  the  network  dynamics  is  in  the  subcritical,  crit¬ 
ical  or  supercritical  domains.  Only  in  the  subcritical  domain  is  the  solution 
to  the  FLE,  that  being  the  MLF,  sufficient  to  describe  the  dynamic  response 
of  an  individual  to  the  other  9,999  dynamic  elements  of  the  network.  In  the 
supercritical  and  critical  domains  the  FLE  is  modified  to  include  a  correlated 
stochastic  driver.  In  the  supercritical  domain  the  relaxation  of  the  fluctuations 
is  exponential  and  solution  to  the  stochastic  FLE  fits  the  individual  autocor¬ 
relation  calculated  data  remarkably  well.  In  the  critical  domain  the  relaxation 
of  the  fluctuations  is  modeled  by  an  inverse  power  law  and  the  solution  to  the 
stochastic  ELE  again  fits  the  individual  autocorrelation  data  surprisingly  well. 

The  lesson  to  be  learned  from  this  combination  of  computation  and  analy¬ 
sis  presented  herein  is  that  complex  networks  of  finite  size  whose  dynamics  are 
members  of  the  Ising  universality  class  described  by  the  DMM  have  an  analytic 
not  just  a  numerical  description.  Instead  of  conhning  the  dynamic  description 
to  that  of  the  macroscopic  variable,  that  being  the  global  or  average  state  of 
the  network,  one  can  also  investigate  how  individual  members  of  the  network 
respond  to  the  influence  of  the  network  as  a  whole.  If  we  consider  the  fluctua¬ 
tions  in  the  global  dynamics  to  be  microscopic,  and  the  potential  of  the  global 
variable  to  be  macroscopic,  then  the  real  time  dynamic  description  of  the  in¬ 
dividuals  is  mesoscopic.  In  general  the  mesoscopic  dynamics  have  a  fractional 
stochastic  differential  representation. 
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